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MARCOS, a numerical tool for the simulation of multiple 
time-dependent non-linear diffusive shock acceleration 
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ABSTRACT 

Wc present a new code aimed at the simulation of diffusive shock acceleration (DSA), 
and discuss various test cases which demonstrate its ability to study DSA in its full 
time-dependent and non-linear developments. We present the numerical methods im- 
plemented, coupling the hydrodynamical evolution of a parallel shock (in one space 
dimension) and the kinetic transport of the cosmic-rays (CR) distribution function 
(in one momentum dimension), as first done by Falle. Following Kang and Jones and 
collaborators, we show how the adaptive mesh refinement technique (AMR) greatly 
helps accommodating the extremely demanding numerical resolution requirements of 
realistic (Bohm-like) CR diffusion coefficients. We also present the parallelization of 
the code, which allows us to run many successive shocks at the cost of a single shock, 
and thus to present the first direct numerical simulations of linear and non-linear 
multiple DSA. a mechanism of interest in various astrophysical environments such as 
superbubbles, galaxy clusters and early cosmological flows. 
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1 INTRODUCTION 

Diffusive Sliock Acceleration (DSA) at supernova remnant 
blast waves is the favoured production meclianism for the 
production of tlie galactic cosmi c-rays (CR) . This theory, 
developed since the late 70s (see lDrurvl[l98^ for a review), 
has now both strong theoretical and observational supports. 
The theo retica l grou nds of the model lie in the early ideas 
of Fermi l| 19491 . 1 1954 ): the regular Fermi acceleration mech- 
anism (known as Fermi I, the stochastic one being known 
as Fermi II) can naturally explain the formation of a power- 
law spectrum by a shock wave - with a remarkable universal 
slope whose value s depends solely on the shock compression 
ratio r (which is always 4 for strong non-relativistic shocks). 
However, the acceleration process can easily be so efficient 
that the CR may back-react on the shock dynamics, modify- 
ing the acceleration process in a fully non- l inear way, and re- 
quiri ng a much more detailed analysis (see lMalkov fc Drurvl 
I2OOII for a review) . 

Thus the DSA mechanism has still received a lot of at- 
tention in the last 20 years, from both a theoretical and a 
numerical perspective. Analytical works have been mostly 
limited to the test-particle (linear) regime. The full non- 
linear time-dependent problem has been mostly investigated 
through num erical simul ations, using several different ap- 
proaches (see |jone3l200ll for a short review). A first class is 



based on particle met hods, from the ear l y Mon te-Carlo sim- 
ulations developed bv lEllison fc Eichler * (1984) to the recent 
Particle-In-Cells codes (eg lDieckmann et al.i,200d ). An alter- 
nate approach consists of solving the (Fokker-Planck) trans- 
port equa tion. This has first been done in the "two fluid" 
model (eg I Jones fc Ka ng 1990l). the n dealing with the com- 
plete particle distribution function (iFalle fc Giddingsl Il987l . 
[Beiilll987l . iKang fc Jon'3ll99ll . iDuffvll 19921 )! 

Most of this work has been aimed at understanding the 
role of single (isolated) supernovae remnants, although in 
many contexts CR are likely to experience many sh ocks, 
most notably inside superbubbles (jParizot et al.l [20041 ). 

In this paper we present a new code for the study of 
DSA, named Marcos for Machine a Accelerer les Rayons 
COSmique^ In section [2] we present the basics of the nu- 
merical methods implemented in our code, which couples the 
hydrodynamical evolution of a fluid with the kinetic trans- 
port of the CR. In section [3] we present the Adaptive Mesh 
Reflnement (AMR) technique which allows us to resolve the 
(very) different scales induced by CR diffusion. In section |4] 
we present parallelization of the code, to be able to study in 
reasonable wall-clock time the effects of multiple shocks. 
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2 ACCELERATING PARTICLES: COUPLING 
HYDRODYNAMIC AND KINETIC 
TRANSPORTS 

In the DSA mechanism particles can be accelerated up to 
very high energies because they are interacting a large num- 
ber of times with a macroscopic structure, namely a shock 
wave (as the forward shock of a supernova remnant). In our 
approach the thermal (fluid being shocked) and non-thermal 
(CR) particles, although intimately coupled, are handled as 
two different populations. The fluid, described by its mo- 
ments p, u, P, obeys the Euler equations, while the CR, de- 
scribed by their distribution function f{x,p), follow a more 
general transport equation. 



2.1 Hydrodynamics 

The hydrodynamics are described by the (non-relativistic) 
Euler equations, which express the conservation of the fluid 
mass density p, momentum pu and total energy density e 
and can be written in the general form: 



dX 
dt 



\7- F X) = 



(1) 



where (from now on in ID slab geometrjjj) the conservative 
variables are 



X 



and their flux 



pu 

,2 



F{X)^ \ pu-' + P 
\ ie + P)u 



(2) 



(3) 



To close this system we consider the usual polytropic 
equation of state of "adiabatic index" 7: 



Pxp'' 

so that the total energy is 
J' ^ 1 2 



7-1 2' 



2.1.1 Hydrodynamic scheme 



(4) 
(5) 



O ur hydrodynamic sche me is adapted from the one used 
in lDownes fc Ravi (|l998f ). It uses an Eulerian finite- volume 
approach, which allows us to take advantage of the conserva- 
tive (hyperbolic) form of the Euler system (see eg iLeVequd 
1 19981 ). In this approach the functions X(x) (where X denotes 



^ For the sake of simplicity the code is presented here in ID slab 
geometry, but it works in ID symmetrical spherical geometry too, 
provided we take into account the correct geometry of the finite- 
volume cells, which are now shells. The formalism of section [2. II 
still holds with u being the radial velocity in equations ([2j and 
^ and with a source term 




any of the three components of X and x is the space vari- 
able) are replaced by piece-wise approximations Xi (where 
i denotes the cell index) . To update the X values from time 
n to time n + 1 we readily discretize equation ^ as 



v.x: 



ViX: 



St 



n+i 



(6) 



where Vi is the volume of cell i, Si_xii {Si+\i2) is the sur- 
face at the left (right) of cell i and F^^^^^j!^ ^® 
flux of X through this surface during the time-step 5t. As 
~ -^{i+i)-i/2 ^'^^^ scheme has the interesting prop- 
erty that (but for boundary effects) the X quantities are 
numerically conservecQ, which is critical for com puting the 
correc t velocity of shock discontinuity. Following iGodunovl 
(| 19591 ) the fluxes F", L are computed using the values of X 
at the interfaces computed by a Riemann solver (we use an 
exact Riemann solver with a linear approximation in smooth 
regions) . 

The code is made second order in time by using a 
leap-frog scheme, and is made second order in space by ex- 
trapolating the X quantities at the cell interfaces before 
running the Riemann solver (doing slope reconstructions, 
which requires "slope limiters" to preserve monotonicity, see 
Ivan Leei)[l977l ). 

We recall that this scheme is subject to the usual 
Courant condition 



^^adv,x ^ 



5x 



(7) 



added to equation ((TJ. 



where timax = {\u\ + c^)^^^^ is the maximum physical signal 
speed at the time considered. 

2.2 Particle Acceleration 

Particles that are energetic enough to be scattered off mag- 
netic irregularities decouple from the thermal plasma (and 
are called CR). It is assumed that this scattering is strong 
enough that the CR are isotropised in both the upstream 
and downstream reference frame. We then describe CR 
through their isotropic distribution function f{x,p,t) de- 
fined so that the CR density is 

n{x,t)= f{x,p,t)4:np'^dp (8) 



and which obeys the following Fokker-Plank equation: 

21 ^ dv£ _ d_ r ^dl\ _l_d^du 
dt dx dx \ dx J 3p^ dp dx 



We note that the Godunov scheme works with the conser- 
vative variables p, pu, e, whereas the Riemann solver works 
with the primitive variables p, u, P. The pressure is computed 
as the difference between the total and kinetic energy of the 
fluid: P = (7 - 1) (e - O.Spn^) so that in strong shocks where 
0.5pn^ S> P and e ~ 0.5pn^ its evaluation might become rather 
imprecise. Actually such a scheme can produce some fake "neg- 
ative pressures" whic h have to be co rrected to some fixed con- 
venient Pm i n (see eg iLeVeguel Il99"3) . To prevent that problem 
iKang et al.l 1 I2OO2I) decided to add another equation specifically 
for the pressure, the equation of conservation of the "modified 
entropy" S = P/p'^~^ (note that this equation is not valid at 
the shock, so that they need to combine it appropriately with the 
usual energy equation). 
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The second l.h.s. term represents convection : CR are ad- 
vected in space as they are bound to the fluid by scattering 
off the magnetic waves present in it. Note that we neglect 
here the movements of the scattering centres (the waves) 
with respect to the fluid (we suppose Va ^ where Va 
is the Alfven speed), so that we don't consider diffusion 
in momentum (known as Fermi II acceleration). The first 
r.h.s. term represents space diffusion of the CR resulting 
from their scattering off the magnetic waves, conveniently 
described by a diffusion coefficient D{x,p) (on which we 
shall focus in section [3}. The second r.h.s term represents 
adiabatic compression of the fluid, the velocity divergence 
^ being the engine of the particles acceleration. 



2.2.1 Kinetic scheme 

From now on we consider that CR are protons (of mass m^), 
and we express CR momenta in rUpC units (and CR velocities 
in c units, CR energy and pressure in r ripC^ units). 

FoUowing iFaUe fc Giddind (|l987h we work with q = 
p^f and y = ln(p) for numerical convenience. We then 
rewrite equation ((9)1 as 



dg_ dug _ d_ ( dg\ _ dg 

dt^ dx ' dx \ dx) "dy ^ " 



d 



dg- 



dg_ 

' dy 



(10) 



where Uy — — appears as a y-advection velocity (note 
also the new source t erm u^q). Our kinetic sc heme follows 
the one presented by iFalle fc Gidd"h^ (|l987l l. but for the 
fact that we use here a Eulerian code so that we have to 
deal with space advection. In fact the particle transport is 
done in two steps using the operator splitting technique. 
Space transport (l.h.s of equation (|9])) is embedded into the 
hydrodynamic Godunov module: CR are transported with 
the fluid as a passive tracer during each hydrodynamic step. 
Then diffusive acceleration (r.h.s. of equation ((9}) is solved 
using a separate finite-difference module. The momentum 
variable y — ln(p) is linearly discretized with step 5y, so 
that in each space cell (indexed by i) we have the full piece- 
wise spectrum of the particles (indexed by j). Ignoring the 
space transport terms, equation (|10|) is discretized as 



5t 



n Trrn . n + 1 i-rr^^ + l 



Sy 



with 



diff" = 



5x^ 



(11) 



(12) 



where we have allowed for a space-dependent diffusion co- 
efficient D to be evaluated at the cells interfaces i ± 1/2, 
and where a;" are coefficients which define the particular 
numerical scheme (see below). 

We first comment on the last two r.h.s terms, repre- 
senting CR energy gain. Y-advection is discretized using an 
upwind scheme, so that j'^ = j and j~ = J — 1 when Uy > 0. 
This leads to a Courant condition 



Sy 



(13) 



dition (|13[) : in that case we just sub-cycle DSA according to 
the hydrodynamic time-step 5^adv,x. 

The Courant condition for the diffusive terms with a 
fully explicit schemes now reads 



StdiS,x < 



2Dn 



(14) 



usually slightly more restrictive than the hydrodynamic con- 



which is much more restrictive than the advection condi- 
tion, because of its quadratic dependence on the space res- 
olution, and because D is an increasing function of p (see 
section I3.1.ip . Hence exploring acceleration to higher max- 
imum momenta requires lowering Stdia,x. To overcome this 
limitation we use implicit schemes which are not limited 
by the Courant conditioi|3 (at the cost of more involved 
computations, and with the risk of loosing accuracy con- 
trol). As seen from equation (lllf) our scheme can be ex- 
plicit (oj" = l,Lu"+^ = 0), implicit {w" = 0,0;"+^ = 1), or 
both, the special case uj" = 1/2, oj"^^ = 1/2 being known 
as the Crank-Nicholson scheme. This scheme is of particu- 
lar interest as it is the only one to be second order in both 
space and time. However it has some drawbacks too, as it 
may give unphysical ne gative values for small values of St 
(|Park fc Petrosianlll996l ). In that respect the fully implicit 
scheme is more robust, but it is less accurate. 

Finally we need to define boundary conditions for space 
diffusion (for space advection CR share the fluid density 
space boundary conditions) : they can be either no flux (re- 
flecting boundaries) or no particle (absorbent boundaries). 
Regarding momentum boundary conditions we simply im- 
pose that p(pmin) = gipmiix) = 0. 



2.2.2 Injection process 

We also need to address the problem of the injection of CR 
from the fluid. As said before the CR are initially noth- 
ing but the high energy particles of the thermal distribu- 
tion. However the problem arising in our two populations 
approach of filling the gap between the fluid and the kinetic 
regimes of the particles is a rather delicate one. 

Injection has usually been parametrized by two quan- 
tities (see eg lFalle fc Giddind Il987l . iKang fc Jonedll99ll ): 
the fraction, rj, of the particles crossing the shock becom- 
ing CR, and the momentum pinj at which they are in- 
jected. The latter can be fixed or defined by pinj = CPth,2 
where pth,2 = y^2m^^fcsTb is the mean downstream ther- 
mal momentum (or alternatively by Vinj = £,'c^,2 where v 
is the particle speed and Cs,2 = ^/ jkBT2/mp is the down- 
stream sound speed) . ^ is expected to be in the range 2 — 4 
(^'/^ — -^2/7 ~ 1.1 for 7 = 5/3). The parameter 77 is less 
constrained, with a typical order of magnitude of 77 = 10""^. 

However iMalkov fc Volkl l|l998t ) have developed a self- 
consistent analytical model of injection of suprathermal 
particles, known as the "thermal leakage" mechanism. 
ICieseler et al.l (|2000l ) have implemented it in the CRASH 
code through the use of a "transparency function" which 
connects the thermal and supra-thermal distributions. They 



* We have also investigated the Supcr-Time-Stepping (STS) 
method which allo ws expUcit schemes to overcome the Courant 
condition (see eg lAIexiades et al.l Il996h . However, this didn't 
prove as time-saving as implicit schemes. 
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have then only one remaining free parameter (the "wave 
amphtude"), which is rather weU constrained - at least for 
strong shocks for which the self-consistent injection rate of- 
ten produces quickly and strongly m odified shocks whic h are 
difficult to handle numerically (see iKang et alj 120021 ') . We 
adopt in our code a simplified re cipe based o n this thermal 
leakage mechanism, proposed bv lBlasi et al.l (|2005l ): given a 
value of ^ we compute rj as 

4 



3^/¥ 



(r-l)C^exp(-eO- 



(15) 



We note that there is also only one remaining parameter, ^, 
which is also well constrained - but rj has a strong depen- 
dence on it. Although very simple to implement, this recipe 
allows a self-consistent description of the leakage of high en- 
ergy thermal particles to the suprathermal population. Note 
that 77 is a function of the shock compression ratio r and 
thus a function of time in the case of a modified shock (see 
section 2. 4p : the r — 1 factor acts as an injection regulator. 
We thus add a source term to the r.h.s of equation 



dF{xs,t) 
dx 



G{x — Xs)5{p — Pinj) (16) 



where a;s is the shock location, F{xs) = {pu)g/mp is the 
particles flux through the shock (evaluated by the code in 

the shock frame), G{x) — -^=-exp^— is a Gaussian 

distribution which spreads the injection around the shock 
location and 5 is the Dirac distribution. Moreover we have 
to take account of the fact that these particles are extracted 
from the thermal population. As usual we neglect the iner- 
tia of the fresh CR, but we remove their energy from the 
fluid: we add a corresponding sink term to the fluid energy 
equation: 



S{x,t)=r,{at)) 



dF{xs,t) 
dx 



G(a 



■ xs 



(17) 



2.2.3 Test-particle acceleration 



In the linear test-particle regime the DSA mechanism is well 
known: theory provides various results we have used to val- 
idate our code and we summarize here. 

Particles injected at constant momentum pinj build a 
power-law spectrum 



f{p) = fop 
of slope 
3r 



1 



(18) 



(19) 



where r is the shock compression ratio (taken as 4 for strong 
shocks). For a strong shock, then, s = 4, for weaker shocks 
s > 4 (that is, / cx is the hardest spectrum achieved 
by linear acceleration by a single shock - without losses) . If 
the injection rate, 77, is constant and the upstream medium 
is homogeneous then the spectrum normalisation is 



fa = s 



■npi 



47rmpp?^j ■ 



(20) 



This spectrum extends from pinj to a maximum momen- 
tum pmax controlled by the scattering of the CR and thus 
in our model by the diffusion coefficient D. If we assume for 
the sake of simplicity that D is constant in space and has 



a simple power-law dependence on p: D{p) — Dop" then p 
grows in time as 



In 



where 

tacc(p) = 



1 , / , t 

- In 1 + a- — ^ 

a \ tacc(Pinj; 



^acc (pi 



Q / 



a = 



- + - ) D{p) 



Ul — U2 \ Ml U2 



(21) 



(22) 



is t he characteri stic acceleration time-scale at momentum p 
fsee lDrurvl[l983l ). 



2.2.4 The non-linear problem: Particle back-reaction and 
shock modification 

It has been noted since the early developements of the DSA 
theory that the CR pressure, defined as 



pv 



/(p)47rp 



47r 



--g{y)dy 



(23) 



grows without limit in the linear regime, which implies that 
some backreaction process occurs. In fact as the CR dif- 
fuse upstream of the shock their pressure gradi ent induces 
a fo rce which pre-accelerates the fluid (see eg iKirkI Il994l 
and iBerezhko fc EUisonI Il999l l : there is formation of a so- 
called "precursor", a smooth and spatially extended ramp 
in the hydrodynamical proflles upstream of the shock. The 
shock itself is thus progressively reduced to a so-called "sub- 
shock", whose compression ratio is Vsub < 4 (however, the 
overall compression ratio rtot from far upstream to far down- 
stream can now be > 4). Particle acceleration is then fully 
coupled with the shock evolution: CR accelerated by the 
shock modify its hydrodynamical structure, modifying in 
turn the acceleration process itself and thus their spectrum 
(which is no longer a perfect power-law, but becomes con- 
cave in shape as particles of different momenta explore dif- 
ferent scales upstream of the shock and thus feel different 
compression ratios). 

Thus we must add the CR pressure contribution in the 
fluid momentum equation as well as in the fluid total en- 
ergy equation (we still neglect the CR inertia, this will be 
discussed at the very end of the article): we add a source 
term 




Q{x) 



on the r.h.s. of equation ([T|. 

Numerically Per is evaluated as 



(24) 



pIj 



3 ^^yrr]? 



(25) 



and its space gradient is computed using centred differences. 
The particle back-reaction is done at each hydrodynamical 
step. If the hydrodynamical and kinetic time-scales are too 
different we may miss the actual CR back-reation. To pre- 
vent the CR decoupling from the fluid we make sure that 
within each hydrodynamical time-step the relative variation 
of the fluid energy Ae/e due to the CR back-reaction is never 
higher than 10%. 



The MARCOS tool for DSA 5 



2.3 Test 1: a strongly modified shocli 

To test our code we start with ti re second test case from 
tlie pioneers of the kinetic approach iFaUe fc Gidding3 (|l987l l 
(hereafter FG). 









"'iH 
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itl 
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40 



2.3.1 Test Design 

The shock wave is generated by a supersonic piston: a pis- 
ton of constant Mach number Mp = Up/cs.i generates a 
shock of constant Mach number Ms = us/cs,i given by (eg 
iLandau fc Lifshit j|l959h 



Ms 



7 + 1 



Mp^ll+l 



7+1 



Mr 



(26) 



We work in the piston frame: the upstream medium is given 
a velocity of —Up, using a "reflecting" left boundary con- 
dition. The piston is then fixed at a; = in the simulation 
box, and the shock emerges out of this left boundary. The 
right boundary condition is chosen to be gradient zero. The 
upstream medium is initially of constant density pi — 1 and 
pressure Pp — 1. The piston Mach number is Mp — 15, so 
that Ms = 100.5, its = 4.52 in the upstream rest frame, and 
r = 3.55. 

There are no CR upstream initially (P°r ~ 0). Par- 
ticles are injected at a constant rate 77 = 0.0225 at a 
variable momentum defined by ^' = 2 (see section r2.2.2|) . 
The piston beta (/3 = Up/c) is adjusted so that initially 
p?jj = 10~^ as done in FG. The momentum grid extends 
from log(pniin) ~ —3 to log(pniax) = +4, with a resolution 
(not critical here) 5y = 0.23 (that is 10 bins per decade). The 
diffusion coefficient is a power-law with a weak momentum 
dependence: D (p) oc p'^'^^. Its normalisation is adjusted so 
that the simulation unit time is the acceleration time-scale 
at injection tacc(pinj) (see equation (|22l) ') as implicitly done 
in FG. 

The simulation is run to tend = 40 as in FG to show con- 
vergence of the coupled fluid-CR system (in the linear case 
we then expect CR to be accelerated to pmax = 3.2 < 4). 
The space box size equals the distance travelled by the shock 
during that time (at constant velocity Ms = 1.77 with re- 
spect to the piston located at x = 0) plus 10 times the 
diffusion length of the highest energy CR (corresponding to 
A = 10 as defined by equation p5p ') that is imax = 250. The 
space resolution is set to Sx = 2, 9.10"^ to achieve numer- 
ical convergence (corresponding to A = 0.050 at pmin (and 
A = 0.037 at p?„j) as defined by equation dMI)- The hydro- 
dynamic Courant number is set to 0.8, the kinetic scheme 
being sub-cycled by another factor of 2. 

2.3.2 Physical results 

The results of the simulation are presented in figures [1] to [l] 
Figure [T] shows the hydrodynamical profiles for various 
times. The shock is visible as a discontinuity traveling to 
the right. As expected (see section r2.2.4|) it is smoothed by 
the presence of a precursor, visible in all profiles, caused by 
the presence of CR upstream of the shock - their pressure, 
added dashed on the bottom plot, exceeds the fluid one. We 
see that after a quick initial adjustment the shock structure 
reaches a quasi-stationary state, as observed by FG. 

Figure [2] shows the time evolution of the shock (in 
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Figure 2. Evolution of the strongly modified shock (in green) of 
figure [T] the shock position Xb, velocity Ms and Mach number Ms 
are plotted versus time. In red are added the results when CR 
back-reaction is turned off, which follow the theoretical evolution 
(dotted). 
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Figure 3. Evolution of some key CR parameters for the modi- 
fied shock (in green) of figure [T] plotted versus time are the maxi- 
mum momentum pcr (its theoretical linear evolution is added dot- 
ted, see equation I I21I I). pressure Per (the fluid pressure is added 
dashed) and adiabatic index 7ci- (the non-relativistic (5/3) and 
ultra-relativistic (4/3) values are added dotted). In red are added 
the results when CR back-reaction is turned off. 
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log(p) 

Figure 4. Time evolution of tlie CR spectrum just downstream 
of the modified sliock of figure [T] Tlie spectrum f{p) is siiown on 
top. On tliis plot the horizontal dotted line marks the theoret- 
ical spectrum normalization /o at injection (in the linear case) 
and the two other dotted lines are power-law spectra of slopes 
s = 4 and si = 4.18 = the theoretical linear slope (for this shock 
compression ratio r = 3.55) and of same normalization /o at in- 
jection. These two remarkable slope values are also marked dotted 
on the bottom plot, which shows the spectrum logarithmic slope 
-s = 91og(/)/aiog(p). 

green). We see that after the quick initial adjustment the 
shock velocity is again almost constant, but slower than in 
the linear case (added in red). 

Figure |3] shows the time evolution of the CR population 
(again with back-reaction switched on/off in green/red). In 
both cases the CR maximum momentum agrees well with 
the dotted theoretical linear prediction (the determination 
of Per is less reliable in the non-linear case, the differences 
seen are not conclusive) . The CR pressure quickly converges 



thanks to the regulation effect of the back-reaction (as both 
the sub-shock mass flux and downstream temperature are re- 
duced) whereas it grows linearly in time in the linear regime. 
The CR adiabatic index is computed as 7cr = 1 + Pci/Ea 
where Per is the " CR fluid" pressure given by equation (|23|l 
and EcT is the "CR fluid" internal energy given by 

£cr = [ K{p)f{p)47Vp^dp ^AH j Xl±Z^g(j/)dy (27) 

J p J y P 

where K{p) is the kinetic energy of a CR of momentum p. As 
expected 7cr starts at the same value as the adiabatic index 
of a non-relativistic fluid 7 = 5/3 (as CR are injected from 
the thermal pool) and goes down (constantly but more and 
more slowly) as CR are accelerated to high energies (tending 
to the adiabatic index 7 — 4/3 of a relativistic fluid). 

The latter quantities are global (macroscopic) proper- 
ties of the CR (seen as a fluid). Figure |3] shows the CR 
spectrum (and its slope) just downstream of the shock at 
some given times. We see how DSA progressively builds the 
CR distribution. Note that the time discretization is lin- 
ear as in figure [TJ so we see that it takes more and more 
time to accelerate particles to higher energies: this is be- 
cause of the diffusion coefflcient dependence on p (see equa- 
tion (|22p ). Initially the CR are injected at log(p['„j) = — 1, 
with a normalization which agrees with the theoretical one 
in the linear case (see equation (|20|l ). But soon afterwards 
Pinj drifts slightly towards lower energies as the downstream 
temperature is reduced (see the bottom plot of figure [l|. 
The CR spectra extend over 20 orders of magnitude, they 
seem to approach a power-law form but the plot of their lo- 
cal slope clearly shows that they are actually concave: they 
are softer at low energies and harder at high energies than 
the theoretical linear slope si = 4.18. This is another well- 
known feature of the acceleration by CR-modified shocks 
(see section I2.2.4|l , due again to the energy dependence of 
the diffusion coefficient 

2.3.3 Comparison with previous study 

The results presented here can be compared to the FG sec- 
ond test (section 4.2, especially their figure 7). Most notably 
the fluid and CR pressure have exactly the same time evo- 
lution, with a convergence at the same values within a few 
percents. The shock Mach number and CR adiabatic index 
also agree well. This is to our knowledge the first direct com- 
parison to the results of FG, and this success cross- validates 
the two codes. It proves that our code can handle well a 
strongly modified shock produced by the tight coupling of 
fluid hydrodynamics and CR diffusive acceleration. 



3 DIFFUSION SCALES AND ADAPTIVE 
MESH REFINEMENT 

We have already seen in the previous section how the energy 
dependence of the diffusion coefficient drives the main fea- 
tures of modified shocks. This is also the reason why realis- 
tic kinetic simulations of DSA are numerically a challenging 
problem: the difficulty is the potentially huge range of space 
and time scales which must be resolved. We now investigate 
this important issue in details, showing its physical grounds 
and its numerical answer. 
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3.1 Diffusion 

3.1.1 The diffusion coefficient 

In our model the scattering of CR off magnetic field fiuctu- 
ations is represented by a diff'usive term (first r.h.s term of 
equation © ) , controlled by the diffusion coefficient D which 
can be expressed as 

D = h,r,{pV (28) 

where v is the particle velocity and /mfp its mean free path 
at this energy. A sp ecial case of interest is t he so called 
"Bohm limit" (see eg lKang fc Jones! Il99ll and lDuffvl[l993 ) 
reached when Zmfp ~ where rg — p/eB is the particle 
gyro-radius, that is when the particles are scattered within 
one gyro-period, meaning that the turbulence causing par- 
ticle scattering is random on the scale rg. This constitutes a 
lower limit on the (parallel) diffusion coefficient, and on the 
acceleration time-scales. This Bohm limit has been widely 
favored in the literature. In that case D oc pv so that 



Db{p) = Do 



(29) 



where one can evaluate Do = 3 x lO'^^ cm^s ^ x 



The dependence of D on p is also frequently conveniently 
modelled by a simple power-law form 



Da (p) = -Dop" . 



(30) 



We note that Db (p) reduces to such a power- law in the non- 
relativistic limit {a = 2) and in the ultra-relativistic limit 
(a = l). 

We can also consider a space-dependent diffusion coef- 
ficient. A common choice is 
Pi 



D{x,p) = ^D(p). 



(31) 



This d ependence on p avoids the sound wave instability stud- 
ied bv lDrurv fc FaUel (|l986h and mimics the compression of 
the magnetic field attached to the density field. 

Note that otherwise the code doesn't handle explicitly 
the evolution in space and time of the magnetic field: some 
constant mean value of B is assumed to get the normaliza- 
tion of D. Therefore we restrict here ourselves to the study 
of parallel shocks. We note that CR are thought to trig- 
ger themselves the field perturbations which make them dif- 
fuse, through various instabilities excited when they stream 
upstream of the shock. Thus it is possible to compute the 
diffusion coefficient D{x,p,t) self-consistently from the CR 
distribution itself, adding the wave transports equations. We 
postpone this problem to a future work. 



3.1.2 The 



fusion scales 



Whatever the precise description used for D, the key point 
is that it is thought to be a growing function of p. As the 
CR momenta span many orders of magnitude (from say p = 
10~^ to p = 10® or even p = 10^ for galactic CR) this 
introduces a wide range of scales. The relevant time scale is 
the diffusive acceleration time-scale (see equation (|22|) 'l: 

D(p) 



The relevant space scale is the diffusion length of the CR 
upstream of the shock: 



t(p) 



Djp) 
us 



(33) 



The space scales range from the microscopic scale where the 
CR decouple from the fiuid (of the order of a few thermal 
gyration lengths) to macroscopic scales (of the order of the 
supernova remnant radius for high energy CR, which then 
escape, thus limiting the Pmax the remnant can achieve). 

From a numerical perspective the resolution of the grid 
is then dictated by the diffusion of the lowest energy CR (we 
must ensure 5x ^ -D(Pmin)/^ts to catch their dynamics well) 
whereas the size of the grid is dictated by the diffusion of the 
highest energy CR (we must SnStir© 2^ max » -D(pmax)/?iS not 
to lose them artificially). The ratio -D(pinax)/D(pmin) (and 
thus the number of cells a:niax/5a;) may exceed ten orders of 
magnitude if D{p) oc p, which is extremely demanding in 
terms of memory requirements and computing time. This is 
the reason why the first simul ations were made with low p 
d ependence of D (a = 0.25 in iFalle fc Giddined (|l987l ) and 
inlKang fc Jonej (|l99ll )). before exploring the Bohm regime 
(|Duffvlll99j 7 - which was achieved by using more involved 
numerical techniques. 



3.2 Adaptive Mesh Refinement 

3.2.1 Principle 

Fortunately we need very high resolution {5x small) only 
around the shock, as this resolution is required by the low- 
est energy particles only, which don't diffuse far away from 
the shock. More generally CR of a given energy require a cer- 
tain space resolution on a certain space extended around the 
shock (the ke y para r neter being a;upst(p)). Hen ce the idea, 
ionee red by buffvl (|l992h and developed by iKang et"al] 
200 ll ). to implement techniques of Adapative Mesh Refine- 
ment (AMR) to allow the numerical resolution Sx to vary 
according to the needs of the CR that are li kely to be found 
at a g iven location at a given time (see also iBerezhko et al.l 
()l994h for a different approach). This allows correct han- 
dling of the transport of CR whereas considerably lowering 
the numerical requirements. 



3.2.2 Design 

We adopt he r e th e technique of nested grids 
(jBerger fc Oligej Il984l l: A'^ sub-grids of increasiiig res- 
olution are added to the base-grid around the shocljj. The 
resolution of the grid at level k (base grid being level 0) is 
Srk ~ Sro/R^ where R is the refinement factor, taken as 
usual to be i? = 2. 

The grid hierarchy is automatically designed by the 
code according to the CR diffusion properties as follows. 
The resolution at the last sub-level A'' is adapted to the low- 
est energy CR, of momentum Pmin = Pn'. 



{5x)j^ = A X a;upst(pjv) 



(34) 



tacc(p) OC 



(32) 



We haven't used tree-based AMR as this technique is much 
more complicated to implement and as its main advantage is its 
versatility but the situation we have to deal with is well defined. 
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where A < 1 (|Kang fc Joned (jigQlh suggest A = 0.05, and 
this is indeed what test 1 of section 12.31 required for full 
convergence) . The resolution of the A'^ — 1 sub-level will nec- 
essarily be 5xn-i — R X Sxn- This resolution will be good 
enough for all particles of momentum above some pjv-i so 
that 5xn-i — X X XupstiPN-i)- Then the sub-grid A'' should 
take care of all CR of momenta from pm = Pmin to pn~i, 
in particular it must contain them as they diffuse from the 
shock, so we set its half-size Ax as 

(Ax)^ = A X XupBt{pN-i) (35) 

with A > 1. The sub-grid A" is now fully defined by its size 
Ax and resolution 5x, the process is then iterated to the level 
A'^ — 1, the resolution of which is imposed by the resolution 
at level A'' and the size of which is imposed by the resolution 
at level N — 2, and so on. The total number of sub-grids is 
adjusted semi-empirically to maximize the AMR efficiency. 

From now on the grids design will be conveniently de- 
scribed by the two parameters A and A defined by rela- 
tions (|34p and 1)35^ . We use the same A for all sub-grids, 
and we extend its definition to a base grid where Ax now 
refers to the distance Xmax — us x t^nd between the position 
of the shock at the end of the simulation and the position of 
the right physical boundary (so that in all cases Ax is the 
minimum distance upstream of the shock). 

3.2.3 Algorithm 

The nested-grids algorithm is a recursive one: at each level 
(from top to bottom) we update the quantities on the whole 
grid, we run the same process at the sub-level, and we re- 
place the grid coarse quantities by the sub-grid finer quan- 
tities (we also correct the coarse fluxes at the sub-grid in- 
terfaces to preserve the scheme conservation properties, see 
iBerger fc Leveau3 Il998l ). We recall that the refinement is 
both a refinement in space (the resolution is divided by R) 
and a refinement in time (because of the Courant condi- 
tion (O). Note that the two operators in the code (the hy- 
drodynamic one and the kinetic one) still operate conjointly 
at each level, as refining them separately would artificially 
decouple their effects. Note also that injection is done at 
bottom level only and propagates to all the upper levels as 
they are updated. 

To set up the refinement the child grid must be given 
appropriate boundary conditions to match its parent pro- 
files. Regarding the hydrodynamics these nested boundary 
conditions consist of a simple filling of the child's ghost with 
the corresponding values of its parent. Regarding the kinetic 
part we note that diffusion is controlled by interface fluxes, 
so that we enforce matching of the diffusion flux between 
the child and parent level. Note that this mechanism works 
whatever the diffusion coefficient scheme (be it explicit or 
implicit), but that the Crank-Nicholson scheme gets more 
sensitive because of such nested boundary conditions. 

The grid hierarchy is set up around the shock position 
at start-up and moves with it over time. For the shock track- 
ing to remain efficient with any number of grids levels we 
allow each sub-grid to move independently both in space and 
time. However a sub-grid can move only by i? = 2 of its cells 
to keep the simple 1 to 2 correspondence of the AMR refine- 
ment scheme (and only at the end of a complete refinement 
step). 



3.3 Test 2: adding Bohm scaling 

Here we extend test 1 using a more realistic (and demanding) 
diffusion model. 

3.3.1 Test Design 

The physical parameters are the same as in 12.3.11 but for 
the fact that we now use the relativistic[f| Bohm scaling 
for the diffusion coefficient, that is D (p) oc p (see sec- 
tion [3TTT}. Numerically we then enter a new world, because 
of the requirements induced by equations (|32|l (longer run) 
and H33p (higher resolution). Simulation of test 1 runs within 

2 hours at high resolution on a desktop-class processor. Now 
with the same Pmax it wouldn't be even possible to allocate 
the grid in memory. We thus apply the AMR technique to 
lower the numerical requirements. We use here A — 0.3 at 
log(pmiii) = —1.5 (that is A = 0.1 at log(pPnj) = —1); and 
A = 6 for each sub-grid (as we have observed that for A ^ 6 
the sub-grids nested boundary conditions for diffusion are 
indifferent for all the CR up to the momentum a sub-grid 
has to deal with) and A = 10 for the base grid. We use a 
better momentum resolution than in test 1: Sy =0.1 (that 
is 23 bins per decade). We've run different simulations with 
different maximum momentum log(pmax) ranging from to 

3 by steps of 0.5. We've run each simulation up to the time 
required for CR to reach this maximum momentum (derived 
from equation (|2ip '). ranging accordingly from 10 to 10000. 
We set the Courant number to 0.5 and we sub-cycle the DSA 
scheme a few times (all the more so since we are at a deep 
level). The number A'^ of sub-grids automatically dumped by 
the code ranges from 1 to 7 as log(pniax) rises from to 3. 

3.3.2 Physical results 

Figures [5] to [8] show the results in the case log(pmax) = 1, 
tend = 100, A'' = 3 (note that on figure [S] one sees how the 
grid hierarchy follows the shock) . The global picture remains 
the same as in test 1: quick convergence to a quasi-steady 
state where the fluid and CR pressures are of the same or- 
der (here the CR pressure doesn't reach the fluid one). The 
effects of the dependence of D on p, which were already vis- 
ible in test 1, are now enhanced: the shock precursor is more 
extended (figure [5]), the CR are more slowly accelerated to- 
wards high energies (figure [S]). 

3.3.3 AMR efficiency 

On figures [S] and [7] the red curves show the results (in the 
case log(pmax) = 1) without activating AMR, that is with a 
single grid having the same size as the base grid (given by 
A — 10) and the same resolution as the deepest grid (given 
by A = 0.3). The results can hardly be distinguished, which 
proves that AMR doesn't compromise the physical accuracy. 
And we show now that on the other hand it considerably 
lowers the numerical cost. Simulations described previously 

^ Using the exact Bohm scaling, that is D (p) oc at injec- 
tion energies, considerably increases the numerical cost of the 
simulation without c hanging much the physical results, see eg 
iKang fc Joiie^ ||2006|) . 
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Figure 6. Evolution of tlie strongly modified shock of figure[5] the 
shock position Xb, velocity Ub and Mach number Mb are plotted 
versus time. The theoretical evolution in the non-modified case 
is added as the dotted line. Results obtained with and without 
using AMR arc shown in green and rod respectively. 
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Figure 5. Time evolution of the hydrodynamical profiles for a 
modified shock (with D{p) tx p, see sections 13. 3. II for simulation 
details). Plotted are the fluid density p, velocity u and pressure 
P (the CR pressure Per is added dashed). The AMR grids hier- 
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Figure 7. Evolution of some key CR parameters for the modi- 
fied shock of figure [5] plotted versus time are the maximum mo- 
mentum Per (its theoretical linear evolution is added dotted, sec 
equation I l21ll 'l. pressure Per (the fluid pressure is added dashed) 
and adiabatic index 7cr (the non-relativistic (5/3) and ultra- 
relativistic (4/3) values are added dotted). Results obtained with 
and without using AMR are shown in green and red respectively. 
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Figure 8. Time evolution of the CR spectrum just downstream 
of the modified shock of figure [5] The spectrum f{p) is shown on 
top. On this plot the horizontal dotted line marks the theoret- 
ical spectrum normalization /o at injection (in the linear case) 
and the two other dotted lines are power-law spectra of slopes 
s = 4 and si = 4.18 = the theoretical linear slope (for this shock 
compression ratio r = 3.55) and of same normalization /o at in- 
jection. These two remarkable slope values arc also marked dotted 
on the bottom plot, which shows the spectrum logarithmic slope 
-s = 91og(/)/91og(p). 



have been made (up to Pmax = 1.5) without using AMR 
too. The computing time of each simulation with (green) 
and without (red) AMR is shown on figure |9] In the case 
log(pniax) ~ 1 presented in section [3.3.21 the speed-up is of 
roughly 250; and it grows steadily as Pmax grows. We note 
that both curves are power laws, with an index 240% times 
lower when AMR is activated. Thus the AMR technique is 
both very efficient and absolutely mandatory to address such 
difficult (realistic) problems. 



12 



Ferrand Downes Marcowith 




1 2 
log (maximum momentum) 



Figure 9. Computing time as a function of the maximum mo- 
mentum pmax (in logarithmic scale) for a Bohm-like diffusion 
(_D(p) (X p). The time unit on this plot is arbitrary, on an Ita- 
nium II processor the simulation to log(pinax) = 1 lasted 13.4 H 
with AMR off and 8 mn with AMR on. The CR are injected 
at log(pinj) = — 1 at t = 0, the minimum momentum is set to 
log(pmin) = —1.5. The two sets of measures have been done 
with (green) and without (red) activating the automatic AMR 
scheme (see sections 13.21 and 13.3.11 for details). Physical results 
for log(pmax) = 1 have been shown on figures [5] to [S] with com- 
parisons of results with and without AMR on figures |6] and [7] 



4 MULTIPLE SHOCKS AND 
PARALLELIZATION 

Even using the numerical trick of AMR, the computing time 
of realistic DSA simulations can still be very high, espe- 
cially if one wants a precise estimate of the CR spectrum 
over a wide range of momenta. This quickly becomes an un- 
acceptable limitation if one wants to investigate the effects 
of multiple, successive shocks. In this section we present an 
original attempt to increase the code computing efficiency, 
aimed at increasing the number of simulated shocks in a 
given reasonable simulation time. 



4.1 Multiple shocks acceleration 

In many astrophysical contexts CR are likely to experi- 
ence many successive shocks: in chaotic stellar winds ijWhitd 
1 19851 ). in rotating accreting flows dSpruit 198^, in radio 
sources with multiple hot spots (Pope et al.lll996h . in the 
galactic center (jMelrose fc Crouch Il997ll. in the OB associa- 
tions inside superbubbles (|Klepach et"alll2000lJ iParizot et al.l 
I2OO4I ). in the early cosm ological flows (^Kang fc Jones| l2005l l 
and in galaxy clusters (|Brunetti fc La zarian 2007|). Many 
efforts have been made to better understand multiple dif- 
fusive shock acceleration (hereafter mDSA), but on quite 
particular cases, and clearly not to the same extent as sin- 
gle diffusive shock acceleration (hereafter sDSA). 



4.1.1 Previous studies 

From a theoretical point of view, mDSA is well understood 
in the linear regime, as we can simply add the effect of a sin- 
gle shock. Being of astrophysical interest it has been investi- 
gated analytica lly since the early developments of the DSA 
theory (see eg lEichleij Il980l . iBlandford fc Ostrikeij Il980l '). 
The main result of mDSA is that the CR spectrum flat- 
tens progressively to a universal asymptotic power-law of 
in dex s = 3 (reg ardless of the shock compression ratios, see 
ee lPope fc Melros e 1994). More involved analytical models 
have been developed including various effects such as second 
order F ermi acceleration, radiation losses (for electrons) an d 
escape (|Schhckeiseij|l984l . lAchterberelll99q . ISchneideilll993l ') . 
But while there is a well-defined analytic framework for the 
linear m DSA, there is n o such thing in the non-linear regime. 
However iBvkovl (|200ll ) has developed a non-linear model of 
the acceleration b y chaotic lar ge scales fluid motions inside 
superbubbles, and iBlasH (|2004 ) has proposed a simple semi- 
analytical model including "seed" particles. 

From a numerical point of view, as we have already 
seen the fully non-linear regime is quite well studied now, 
but in the single shock model: although feasible, multi- 
ple DSA has rec e ived e xtremely reduced attention so far. 
ICieseler fc Jonesl (|2000l) have studied acceleration at mul- 
tiple oblique shocks, but in the test-particle regime (they 
found the CR spectr um to harden substantially). Recently 
iKang fc Jone3 (|200i) have investigated the effect of an up- 
stream CR pressure in the non-linear regime, but for single 
shocks (they concluded that it doesn't affect strong shocks 
much but may enhance the efficiency of weak ones). 

Our aim with the code presented in this paper is to be 
able to study CR acceleration by multiple shocks in detail, 
in the full time-dependent non-linear regime. 



4.1.2 Inter-shock physics 

A new important point to consider when simulating multiple 
shocks is the fate of the CR between two successive shocks. 

First the shocked fluid will decompress to recover 
its initial state, and the CR being bound to it through 
scattering will exper ience adiabatic decompression (see 
iMelrose fc Popel Il993l ) : when the shocked fluid density is 
decreased by a factor r the CR momenta are decreased by 
a factor r^''^ . We would like to insist on the physical impor- 
tance of this decompression, whose inclusion is essential in 
the correct treatment of mDSA. For instance a sequence of 
identical shocks of ratio r will produce in the linear regime a 
CR spectrum with the well-known s = 3 slope if and only if 
the CR are decompressed by the corresponding r^^^ factor 
between each shock: if we don't decompress them enough 
they will pile-up from the injection momentum, producing 
much harder spectra, and if we decompress them too much 
they will still form a power-law but a steeper one of slope 
s > 3. 

Apart from these energy losses due to decompression 
the CR might simply escape the system before the next 
shock occurs. In fact CR can be quite well conflned in a 
medium thanks to their diffusion, whose typical length scale 
grows with the CR energy, so that the CR spectrum will be 
depleted from its highest part. In a given physical situation, 
given the typical time between two shocks and the diffu- 
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sion experienced by CR during that time (not necessarily of 
Bohm type) we can estimate the maximum momentum of 
the remaining CR when the next shock arrives. 

We also note that the second order Fermi acceleration 
mechanism (diffusion in momentum) might become an im- 
portant process between the shocks However we won't in- 
vestigate this possibility any further in this paper. 



4.1.3 Numerical treatment 

The code is fully automated to run a sequence of shocks. 
Assuming that the fluid has enough time to decompress be- 
tween two shocks, the same initial hydrodynamical condi- 
tions are used for each shock. The shock Mach number (and 
thus velocity and compression ratio) can however vary from 
shock to shock. At the end of each shock we take the down- 
stream CR spectrum, modify it to take into account the 
inter-shock physics of section l4.1.2l and pre-inject it in each 
space cell before launching the next shock. Although this 
mechanism seems simple we have to elucidate a few points. 

First the very idea of a "downstream spectrum" sup- 
poses that we have reached so me converged state d own- 
stream of the shock. According to lKang fc Jonej (|2005l ) such 
a quasi-stationnary state is obtained once the CR are accel- 
erated up to the relativistic regime (p > vripc), which agrees 
with our own observations. In the following we consider that 
the time each shock is run is long enough so that each sin- 
gle shock can fully relax regarding particle acceleration (at 
least up to the maximum momentum we consider), so that 
the downstream CR pressure is well defined. 

To mimic adiabatic decompression the CR spectrum 
ln(/) is shifted in y = ln(p) by 

Ay(r) = i ln(r) (36) 

towards lower energies. Lost values of / below ymin are sim- 
ply discarded. We have checked by lowering the value of 
2/min that they don't influence the overall subsequent spec- 
trum evolution. Missing values of / between i/max — Ay and 
J/max are filled by linearly extrapolating the slope from the 
end of /. This treatment of the gap gives the best accu- 
racy at high energies in the linear regime. As emphasised 
in section [4.1. 2l to obtain correct results we need to resolve 
precisely this decompression shift and thus to use a high 
resolution in momentum: 5y <^ Ay. We have found that in 
order to obtain exactly s = 3 in the linear regime, the shift 
Ay/Sy must in fact be an exact number of bins, and must 
be as high as roughly 10. But Ay depends on r, and we 
want to be able to run multiple simulations with variable r, 
and in non-linear simulations r will be constantly modified 
(see section [2. 2. 4|) so that we will never know its final value 
beforehand. To solve this problem we proceed as follows. 
We use the same momentum resolution Sy to run all shocks, 
fixed so that 5y < Ay{rmin) / J where J is a chosen integer 
^ 1 and rmin is the lowest allowed compression ratio (note 
that in the case of modified shocks the relevant ratio for de- 
compression is the total one, which will always be greater 
than the ratio imposed initially, so that rmin is well defined) . 
At the end of each shock we measure the actual value r of 
the compression ratio, re-bin the numerical spectrum / with 
new resolution Sy' = Ay{r)/J, shift this under-sampled / 



by Ay(r) that is by exactly J bins, and then re-bin / back 
to the nominal resolution Sy. 

Regarding escapes due to losses we simply give to the 
code a cut-off momentum pcut above which the CR spectrum 
/ is set to zero. 

4.2 Parallelization 

Even using AMR running multiple realistic shocks simula- 
tions can easily be very demanding in computing time, lim- 
iting drastically the possibility to do parameter studies and 
thus fully explore the mDSA mechanism. We note here that 
even if the code is ID in space the inclusion of the full spec- 
trum of particles in each cell makes it actually 2D. Thus we 
now pay more attention to the momentum dimension. 

4.2.1 Principle 

Confronted by this problem, Ijones fc Kangj (|2005h use an 
interesting " coarse-grained finite momentum volumes" tech- 
nique to lower the constraints imposed by the p dimen- 
sio n. The basic i dea u nder t his approach (first introduced 
by lJun fc Jone^ (|l999l ) and Ijones et all (|l999l )) is simply 
to lower the numerical resolution Sy in momentum, but pre- 
scribing a power-law spectrum shape to each part of the dis- 
cretized spectrum in order to keep reasonable accuracy. The 
numerical spectrum is then no longer a piece-wise constant 
function but a piece- wise linear function. This technique al- 
lows reasonably good estimates of the modified shock evo- 
lution with unusually low momentum resolutions (which 
can be as low as 2-3 bins per decade). However we believe 
that the adiabatic decompression between multiple shocks 
wouldn't be handled properly by such low resolutions, as it 
is typically of only 1/5 of decade and has to be well sampled 
to get precise results in the linear regime (see section l4.1.3p . 

In this paper we adopt an alternative way to lower the 
p— dimension numerical cost, without any compromise re- 
garding momentum resolution, which consists simply of fully 
exploiting the power of modern super-computers by paral- 
lelizing the code. We consider here the paradigm, imple- 
mented with MPI, which consists of splitting the grid over 
many processors, each processor running the same code but 
on its own data - which still involves some communications 
between the processors, most notably to define their bound- 
ary conditions. 

4.2.2 Implementation 

We have parallelized our code in the momentum space as 
this is straightforward to implement (as nothing happens in 
momentum space but a global advection to higher energies) 
and perfect load-balances can always be obtained (provided 
we slightly adjust the resolution Sy so that the number of 
momentum bins is an exact multiple of the number of proces- 
sors). Parallelizing in space would be more difficult because 
of space diffusion, and far less efficient because with AMR 
realistic problems are always ill-balanced, no matter how we 
split the grids. However parallelization in momentum suffers 
from two limitations. First not all the code is parallelized 
but only the parts dealing with CR and the maximal ef- 
ficiency of the parallelization of a code is always limited 
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by its sequential portions - but realistic problems are CR- 
dominated so that it is easy to reach very high ratios of the 
parallelized over un-parallelized fractions of the code. Sec- 
ondly, this ratio determines only the maximum acceleration 
achievable through parallelization: in practice the effective 
efficiency of parallelization is limited by the extra cost in- 
duced by inter-processor communications, which eventually 
leads to a saturation of the efficiency - but we managed to 
achieve good scalings as shown below. 



4.3 Test 3: doing multiple shocks 

Here we present the evolution of test 2 when multiple shocks 
are run. To the best of our knowledge, these are the first 
direct simulations of time-dependent linear and non-linear 
multiple DSA. 



4.3.1 Test Design 

We start from the design of section 13.3.11 with the same 
hydrodynamical initial conditions for each shock, but an 
evolving CR population, as explained in section 14.1.31 We 
recall that particles are injected at log(pinj) = —1. We want 
here to study acceleration up to log(p) — 1. When multiple 
shocks are run with CR now present everywhere upstream 
the code faces harder numerical precision issues at high en- 
ergies. To fix that first we set a bigger maximum momen- 
tum log(pinax) = 1.5 and we set bigger sub-grids (A — 15). 
We use the same space resolution as before: A — 0.3 at 
log(pniin) ~ —1.5 (that is A = 0.1 at log(p?jj) = —1). We use 
a better momentum resolution; Sy — 0.036 corresponding 
to exactly 64 bins per decade (which, over 3 decades in p, 
allows perfect load balancing when running the code in par- 
allel on clusters of processors that are powers of 2 up to 64) . 
This nominal resolution is adjusted at the end of each shock 
to be exactly 10 bins per decompression shift (as explained 
in section I4.1.3|l . We consider here that pcut > Pmax so that 
CR don't escape between two shocks ( section l4.1.3p . 



4.3.2 Parallelization efficiency 

Figure [TO] shows the gain brought by parallelization. We ob- 
tain good scalings up to a few tens of processors. Thus par- 
allelization allows us to study the acceleration by multiple 
shocks within the time previously required to study acceler- 
ation by a single shock. 

In this simulations the DSA operator represents a bit 
more than 90% of the computations time, and as the hy- 
drodynamic operator also advects CR the fraction of the 
code actually benefiting from parallelization is more than 
99.5%. On slightly less CR-dominated simulations we have 
observed that the scaling is better on the shared memory 
machine than on the distributed memory machine, as our 
code is then bound by communications. The slightly less 
good scaling observed with 64 processors is no surprise given 
that in that case each processor deals with only 3 momenta 
cells, which makes a high surface/ volume (that is communi- 
cations/computations) ratio. We note here that a good point 
of parallelization in momentum is that it's all the more use- 
full since one wants to investigate high energies. Indeed the 




number of processors 



Figure 10. User computing time as a function of the number of 
processors for simulations of section 14.3.11 Measures have been 
made on two machines, the soleil super-computer of the French 
Calmip collaboration (a cluster of 120 Itanium II processors with 
shared memory) and the rowan super-computer of the Irish Cos- 
mogrid (a cluster of 256 Xeon processors with Gigabit Ethernet). 
Tests have been made up to the maximum number of processors 
available on this middle-class machines (32 on soleil and 64 on 
rowan). The computing times have been normalized to empha- 
size the parallelization scaling with the number of processors. The 
three dotted lines show the theoretical scaling for a perfect paral- 
lelization (that is with no induced over-cost) of respectively 99% 
and 100% of the code. 



higher Pmax, the bigger the momentum grid, and the more 
processors one can use with a same given efficiency. 



4.3.3 Linear results 

Figures [11] and [12] show the evolution of test 2 when 30 
such shocks are successively launched with CR back-reaction 
turned off. In this linear case the CR pressure grows forever, 
so that we end each shock at the time tend = 300 correspond- 
ing to log(pniax) = 1.5. Figure [TT] shows the evolution of the 
CR spectra just downstream of the last shock. The spectra 
are all normalized at the injection momentum to emphasize 
the slope evolution (the actual normalization rises by a fac- 
tor of roughly two). We clearly see the convergence of the 
spectrum from an initial power-law of slope s — 3r/{r — 1) 
(the well-known linear solution for a single shock) to a final 
power-law of slope s = 3 (the well-known limit in the case 
of multiple shocks). The slopes plot (bottom) shows that in 
between the spectrum is never a simple power-law, as the 
asymptotic convergence to s = 3 is all the more slow since 
the momentum is high. The way the spectrum hardens at 
different momenta is shown on figure [T^ where we plot the 
slope as a function of the number A'' of shocks. Dashed is 
added the theoretical slope computed (with centered finite 
differences) from the analytical expression of the spectrum 
produ ced after N shocks, which reads (eg iMelrose fc Popel 

mi)) 
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Figure 11. Time evolution of the final downstream CR spec- 
tra for a sequence of successive linear shocks (see sections 14.3.11 
and l4.3.3l for details). Each coloured line shows the CR distribu- 
tion just at the end of a shock. The spectra /(p) are shown on top 
(all normalized so that /(pinj) = l)i where we have added (dot- 
ted) the three power-laws of slope s = 3 and s = 4 and si = 4.18 
the theoretical linear slope for the compression ratio r = 3.55 of 
the shocks. These three remarkable slope values are also marked 
(dotted) on the bottom plot, which shows the spectra logarithmic 
slopes —s = 9 log(/)/9 log(p). The evolution of the slope from si 
to 3 with the number of shocks is shown on figure 1121 for three 
different momenta. 
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where R = 



r.l/3 



is the decompression factor and si = 



3r/(r — 1) is the single shock slope. The code reproduces the 
expected behavior within roughly 1% for the three momenta 
(close to the maximum momentum the code is less precise). 
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Figure 12. Evolution of the final CR spectrum slope at three 
different momenta (log(p) = —0.5, 0.0, -1-0.5) for a sequence of 
successive linear shocks (this figure consists in three vertical cuts 
in the bottom plot of figure llll l. The theoretical results (from 
equation i37l ) are added dashed. The three remarkable slopes 
s = 3, s = A and si are marked dotted. 



This validation of our code in the linear regime gives us 
confidence to explore the unknown non-linear regime. 



4-. 3-4 Non-linear results 

On figure [13] we now allow the CR to back-react on the 
shock. Each shock is run until the downstream CR pressure 
has converged, before doing decompression and launching 
the next one. We recall that one of the consequences of CR 
back-reaction is that pinj varies in time (as the downstream 
state changes). However we still use here a fixed injection 
fraction rj, but we run different simulations with different 
values of as this parameter is quite poorly constrained. To 
get an overall picture we use a broad range of values of rj, 
around the value 770 = 0.0225 used since test 1 to match FG 
parameters. For 77 ^ 10~^ (a value which seems physically 
unreasonably high) the very first shock gets fully smoothed 
by CR back-reaction before Per has converged (following FG 
we consider a shock to be smoothed - and thus stop injection 
of fresh CeQ - as soon as the Mach number of the sub-shock 
drops below Ms, cut ~ 1.3). As 77 is lowered to around 10^^ 
the first shock can run but the cumulative effect of shocks 
is such that one shock eventually gets fully smoothed. As 
77 is lowered to around 10~^ the number of shocks A'" be- 
fore smoothing occurs raises exponentially: below roughly 
77c — 1.5 X 10~^ it seems that virtually any number of shocks 
could run (although all this shocks are still modified ones). 
We have limited here the maximum number of shocks to 30, 



^ Note that stopping injection doesn't mean stopping accelera- 
tion. FG have shown that it is possible to have a CR-dominated 
"self-sustaining" shock without injection of fresh particles. They 
had also no advection of particles, whereas with multiple shocks 
we have to consider the effect of an upstream population too. We 
postpone the detailed study of this aspect to a future work. 
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as this seems reasonable for both numerical reasons (given 
the evolution of the red curve it would take extremely long 
times to fully explore the very low injection fractions), and 
physical reasons (considering a few tens of successive strong 
shocks makes sense in environments such as superbubbles). 
At the last "complete" shock N we measure the range of 
spectra slopes s (between the injection momentum and the 
momentum of hardest slope, at which its final decay starts). 
We observe two evolutions as rj decreases. First the spectra 
globally harden as rj is lowered, which is expected as more 
shocks can run. Note that as we limit ourselves to A'' = 30 
the slopes below rjc can't be directly compared with the 
slopes above rjc'. below rjc the slopes would get closer to the 
s = 3 limit if one would allow for a higher number of shocks. 
Anyway we see that in the non-linear regime the building of 
the s = 3 spectrum within 30 shocks (as on figures [Tl1 and[T2)) 
requires an injection fraction lower than ri = 10"''. Second 
the range of slopes gets constantly narrower, especially be- 
low rjc (that is when when CR no longer limit the number of 
shocks). Thus this simulations suggest the existence of two 
regimes of mDSA with respect to the injection ratio 77: there 
seems to be some critical rjc (here of roughly 1.5 x 10"^) 
above which CR dictate the fate of the shocks (producing 
soft and irregular spectra) and below which CR are almost 
transparent to the successive shocks (producing harder and 
more regular spectra). We have observed the same global 
picture with other simulations (not shown here) involving a 
constant diffusion coefficient D. 

We note that the self-c onsistent injection fraction pro- 
posed bv lBlasi et al.l l)2005i 'l (equation (|15p l is here initially 
rjB ~ 10""*^, thus in the regime where CR dominate from the 
very first shock. This self-consistent rjB is time-dependent 
and is lowered as the shock gets modified, but we have ob- 
served that the first shock still gets fully smoothed before a 
quasi-steady state has been reached. Such a very high back- 
reaction might be surprising for a thermal leakage mech- 
anism. It comes from our particular choice (to match FG 
parameters) of the ratio of the velocity of injected CR to 
the downstream sound speed ^' — 2, as rjB has a very strong 
dependence on this free parameter (recall that ^' = 1.1^). 
^' = 2 is a realistic but rather low value, we could suggest as 
well — 3, in which case rjB is initially ~ 270 times lower, 
that is rjB ~ 3.7 x 10~*, that is in the regime where CR 
are transparent to the successive shocks. Thus this points 
out that Blasi's model, although providing a self-consistent 
injection fraction, still requires some initial tuning. 

Finally we note that the actual fate of CR-dominated 
shocks might depend on geometry effects too. For the sake 
of (both numerical and physical) simplicity we have consid- 
ered here piston-driven shocks in slab geometry, this work 
shall now be extended to supernova-like shocks in spherical 
symmetrical geometry. We also recall that all this results 
were obtained with a numerical momentum box limited to 
log(pinax) = 1.5. Thus their validity depends the assump- 
tion that CR accelerated to higher energies: (i) don't have a 
major impact on the final shock structure (and thus on the 
spectrum shape at lower energies); (ii) have enough time to 
escape from the system between two shocks. Another more 
fundamental assumption on which we rely (as all other stud- 
ies of that kind) is that the inertia of the CR is negligible. 
Although reasonable for sDSA this might be questionable for 
mDSA because of the cumulative effect of shocks. However 
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Figure 13. Evolution of the number of shocks A*^ (up to 30) that 
reach a quasi-steady state before complete smoothing and of the 
range of final CR spectra slopes —s as a function of the injection 
fraction jj in the non-linear case (see section [4.3.41 for details). 



we note that the higher the injection fraction, the lower the 
number of shocks that we actually run. And we have checked 
that the ratio Pa/p remains always below 10 % in all our 
non-linear simulations. 



5 CONCLUSION 

We have presented a new code aimed at the simulation of 
time-dependent non-linear diffusive shock acceleration. It is 
based on the kinetic approach, coupling the hydrodynamical 
evolution of the plasma with the diffusive transport of the 
distribution function of the supra-thermal particles. As such 
it falls under the legacy of the pioneers (Fallc & Giddin, 



19871. lDuffvlll99d) and of the masters (Kang fc Jones 19911 
Kang et allfeoOll ) of the genre. As the CRASH code it im- 
plements an efficient AMR technique to deal with the huge 
range of space- and time-scales induced by CR diffusion of 
Bohm-like type. To save even more on computing time we 
have also parallelized our code in momentum so that we can 
study acceleration by multiple shocks as fast as accelera- 
tion by a single shock. However in many aspects (high-Mach 
flows, shock tracking, self-consistent injection) our code re- 
mains numerically simpler than CRASH - which can be 
both a limitation and an advantage. Regarding the physics 
we note that various mechanisms of importance could be 
included in the code: self-consistant diffusion coefficient 
(adding magnetic waves transport), second-order Fermi ac- 
celeration (especially between multiple shocks), electrons ac- 
celeration (adding radiative losses), CR radiation (hadronic 
and leptonic). . . 

We have presented a few tests that show that our code 
works well, with respect to both the physical accuracy and 
the numerical efficiency, even in realistic difficult situations. 
We are now able to investigate in details the various as- 
pects of the DSA mechanism, which 30 years after its early 
developments still poses some difficulties. In particular we 
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can address the non-linear multiple DSA mechanism, which 
we believe hasn't received so far all the attention it deserves. 
Our very first results suggest that the injection fraction plays 
a crucial role. We intend now to study in more details the 
situation in superbubbles. 
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